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ABSTRACT 

Thermal surface radiation has been detected with X-ray instruments for several neutron stars 
with high spectral, spatial and timing resolution. These observations allow for direct study of 
fundamental properties of the source, but require model atmospheres and spectra for careful 
interpretation. We describe a robust and extensible implementation of complete linearization 
for computing the spectra of isolated cooling neutron stars for a broad range of temperature 
and magnetic field. Self-consistent spectra are derived for arbitrary magnetic field geometries 
at B < 10 14 G. 
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polarization - magnetic fields 



1 INTRODUCTION 



Isolated neutron stars (NSs) are most commonly discovered as ra- 
dio pulsars and only infrequently as soft X-ray sources or through 
associations with supernovae. Cooling of the hot, isothermal NS 
^ ■ core (T c ~ 10 9 K) proceeds initially by neutrino diffusion, fol- 



lowed by an extended period of thermal emission at soft X-ray en- 
ergies. In this latter phase, conduction in the insulating iron crust 
gives way to photon diffusion and free propagation of radiation 
through the tenuous outer layers of the NS atmosphere. The thermal 
evolution of NS has been the subject of much theoretical considera- 
tion si nce the seminal work of lTsuruti]! 19641) . and lChiu & Sa lpeter 
ll964l) . In principal, measurement and interpretation of thermal X- 
rays emitted at the NS surface provides direct evidence of the star's 
surface composition, magnetic field properties, and, in the absence 
of reheating, may constrain the cooling history of the NS popula- 
tion or the equation of state of the stellar core. In practice, detection 
of these thermal X-rays is difficult. The stellar environment is en- 
ergetic, particularly for young NS, and non-thermal emission pro- 
cesses which originate in the active pulsar magnetosphere dominate 
the surface radiation. Young pulsars may be obscured by diffuse 
non-thermal X-ray emission from the supernova remnant or from 
bright compact nebulae (plerions) driven by pulsar winds. Some 
NS are exceptionally luminous in thermal X-rays possibly owing 
to either a highly transparent surface composition or to reheating 
processes which maintain their core temperature, while the thermal 
and (brighter) non-thermal emissions for some cooling NS decline 
presumably at different rates until the surface radiation can be dis- 
tinguished at late times. 

The fundamental properties of NS thermal radiation are 
regulated by the source luminosity, the strong surface grav- 
ity of the star, and the composition of the atmosphere which 



is likely determined by either accretion (from the interstellar 
medium, a fossil debris disk or stellar companion) or fall-back 
of supernova ejecta. Provided accretion is neither too great nor 
continuous, gravitational settling of heav y elemen ts occurs on 
short timescales iBrown. Bildsten & Rutledgelll998h . leaving the 
lightest elements to form the X-ray photosphere; only a tiny 
fraction of the stellar mass in these light elements is neces- 
sary to form an optically thick layer. Spectral reprocessing by 
the light element plasma atmosphere produces surface emis- 
sion which differs substanti ally from the blackbody function re- 
gardless of magnetic field iRomani 1987; Raiasopal & Romani 
| l99dlZavlin. Pavlov & Shibano\il996bHLlovd. Hernauist & Hevll 
2003), and these discrepancies motivate the accurate computa- 
tion of model atmospheres and spectra. Tabulations of unmagne- 
tized NS X-ray spectra are now available in the XSPE C soft ware 
IPavlov et alJll992HZavlin et alJll996bl : lGansicke et alj|2002t) . al- 
lowing for immediate comparison with the blackbody model for 
B < 10 10 G. The optical and near UV thermal flux shows strong 
dependence o n both composition IPavlov et alll996 ) and magnetic 
field strength iLlovd et all2003l) . The most persuasive models for 
NS thermal radiation must reconcile the X-ray spectrum with opti- 
cal and UV fluxes whenever possible. 

Several codes have been developed to calculate the atmo- 
spheric structure and associated thermal spectrum for a variety of 
NS surface compositions and magnetic field regimes. The atmo- 
sphere is described by the self-consistent solution to a system of 
radiative transfer equations and equilibrium constraints. Essentially 
all methods for calculating radiative transport in magnetized NS at- 
mospheres proceed from the work o f Gnedi n & Paylovl < 19751) . who 
established the connection between the density-matrix formalism 
and the (more immediately tractable) formulation of radiative trans- 
port in two coupled "normal" polarization modes; i.e., the limit of 
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large Faraday depolarization. The plasma is strongly polarizing to 
X-ray radiation at magnetic field strengths comparable to those of 
radio pulsars. The optical and UV spectra of such NS are also po- 
larized for modest B, but measurement at these wavelengths is dif- 
ficult owing to the presence of non-thermal radiation, absorption 
in the interstellar medium, and the overall faintness of the sources. 
The large plasma anisotropies induced by magnetic fields strongly 
affect the angular distribution of the emergent radiation field of the 
atmosphere, and a complete model of the thermal X-ray emission 
of a typical NS requires integration over intensity profiles from sur- 
face elements having local temperatu re and magnet ic field which 
vary across the stellar surface (Zav lin et alJll995al) . All modern 
strategies appeal to the full angle and energy dependence of the 
plasma opacity and coupled radiation field. 

Early efforts to model the X-ray emission from NS atmo- 
spheric plasma were preceded by the development of the radiative 
transfer by the nor mal mod e approximation in uniform (isother- 
mal) plas ma slabs {Kam inker et all l 19831) using the normal mode 
strategy of Gnedin & Pavlov] 1 19751) to describe the radiative dif- 
fusion. IShibanoy et alj 119921) applied the diffusion method of 
iKaminker et alJll983h to NS atmospheres with fully ionized H/He 
mixtures by implementing the constr aints of hydrostatic and radia- 
tive equilibrium. IPavlov et alJ 1 19941) computed spectra for gener- 
ally oblique magnetic fields the NS atmosphere from the thermal 
structures from lShibanov et aT] J 19921) by recalculating the plasma 
opacities and making a formal integration of the transfer equation 
from the prior solution. IPavlov et al J JT995) improved the diffusion 
results by using them as an approximate s olution to the temperature 
correction procedure o f A tier & Mihalasl Jl968l) : this work also in- 
cluded a prescription for th e ionization equ ilibrium of the Hydro- 
gen plasma (see also lZavlin & Paviovll2002l) ). 

iRomanil ll987l) described the first self-consistent NS model 
atmospheres for weak magnetic fields using the Los Alamos Opac- 
ity Library and some simplifying assumptions about the relative 
contribution of thermal and scattering emissivities. In this method, 
the radiation field is derived from application of the Milne in- 
tegral operator to the source function and temperature correc- 
tions to the atmospheric struct u re were ev aluated by the L ucy- 
Unsold method iMihalas 1978). Raiagopa r& Romanil 1 19961) re- 
visite d this work using the O PAL tabulations llelesias et all 19921) 
while iRaiagopal et alJ 1 199% . using essentially the same formal- 
ism, extended their opacity model to roughly approximate the 
bound state transitions in magnetized iron, and included electron 
scattering and free-free absorption in the normal mode approxi- 
mation. Each of these efforts revealed the extent of compositional 
effects on the thermal X-ray spectra of H e, Fe and cosm ic abun- 
dance atmospheres. More recently, Gansicke et al. 120021) have re- 
visited the low-field regime using these methods to consider reheat- 
ing in the outer atmospheric layers of some NS. IWerner & Deetienl 
(2000) have taken a different computational approach, deriving the 
first non-LTE model calculations in the weak field regime using the 
Accelerated Lambda Iteration combined with the Opacity Project 
cross sections; these authors find modest corrections due to NLTE 
effects in light element plasma at low temperatures, but these are 
potentially more significant for magnetic fields stronger than they 
consider. Some recent modeling efforts have specialized to the cal- 
culation of light element "magnetar" spectra (B > 10 14 G) which 
include some prescription for vacuum polarizat ion effects which 
are particularly relevant for B > 4.4 x 10 13 G. IZane et aTlfeOOOl) 
considered in particular atmospheres heated externally at low ac- 
cretion rates, generalizing energy balance in the atmosphere to 
account for energy deposition by charged particles. \Ozej 1200 ll) 



and H o & Lail 1200 ll) have each generalized the the Lucy-Unsold 
method to the normal mode approximation to derive model spectra 
for highly magnetized, light element NS plasma. Ihq & Lail 1 2001 ) 
also considered the diffusion equation solutions to perpendicular 
(normal) and parallel magnetic field geometries in the intense field 
regime. 

This article describes an implementation of the classical tech- 
nique of complete linearization (CL) to the modeling of elementary 
NS atmospheres in radiative equilibrium. We restrict our derivation 
to model atmospheres with pure Hydrogen composition in the limit 
of complete ionization, but the formulation is immediately useful 
to admixtures of Helium. The method is robust, flexible with re- 
spect to the range of model parameters for which it can be applied, 
and is extensible to detailed statistical equilibrium for partially ion- 
ized models and complete Stokes transport. Self-consistent model 
results for magnetic fields B ^ 10 14 G and of arbitrary field ori- 
entation, and a range of flux temperatures characteristic of NS ob- 
served in thermal X-rays have been calculated with this method- 
ology. The article is organized as follows: The magnetized plasma 
opacity is derived in [J2] with attention to the role of protons and 
vacuum effects in radiative processes. The radiative transfer equa- 
tions and their constraints are derived in J3| The grid methods and 
solution techniques are described in An abbreviated selection 
of model results are presented in Sj5| and emphasize the angular 
intensity properties and spectral characteristics, including polariza- 
tion for the weak, intermediate and strong field regimes. Extensions 
and limitations of the present method are outlined in J6| 



2 OPACITY AND EMISSIVITY OF MAGNETIZED 
PLASMA 

The detailed polarizing properties of a light element magnetized 
plasma can be derived from its dielectric response to electromag- 
netic w ave s. This presentation i s simil ar to that found in Ventura 
1 1979) and Mesz aros & Venturi] 1 19791) . The equation for an elec- 
tromagnetic wave propagating in a medium of refractive index N 
along the direction k, including the vacuum permeability, is 



-aN 2 (l- kk) + hN 2 {k x b){k x b) + e E = 



(1) 



where b = (sin 6, 0, cos 9) is the unit vector in the direction of 
B, and a and h are defined below. For uniformly magnetized cold 
plasma, the vacuum corrected dielectric in the B || z frame is 
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in which the plasma and cyclotron frequencies for each species en- 
ter through the parameters 
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for u>p 3 = 4ne 2 n B /m s and lj Cj3 = q s B /(m s c); note that vhj 2 
is a signed quantity. The damping term A s = 1 + iu s /u splits the 
plasma dielectric into hermitian and anti-hermitian contributions, 
but when evaluating the polarization of propagating waves it is ac- 
ceptable to consider only the hermitian part of the dielectric. The 
parameters a, h and q are functions solely of the magnetic field 
strength: 



a(x) = 1 + — [-2X (ar) + xXfa)] 
h ( x ) = ^ [x 2 Xq(x) - xXq(x)] 
q(x) = 



(8) 
(9) 
(10) 



where x = Bq/B, Bq 



'/(eh) 



4.4 x 10 13 G is the 



quantum critical field strength, and the functions Xq and Xi are 
found in iHevl & Hernauistl 1 19971) : the expressions J81 1 01 are valid 
for all field regimes and revert to the results of Meszaro s & Ventur3 
in the weak field limit, x > 1. The tensor e in the k || 
z frame is recovered from J5J by an orthogonal transformation of 
angle cos 9 = k ■ b about the k x b axis. 

The polarization vector for an EM wave in mode j in rotating 
coordinates about the magnetic field is 

e J 1 = ^-i= [K 3 X cos 9 - K{ sin 9±i], K 3 Z cos + K j x sin 6^ (11) 

where the component indices are (±, 0), respectively. The normal 
modes are usually referred to as ordinary (O) and extraordinary 
(X); waves propagating in the ordinary mode are largely unaffected 
by the magnetic field, while X mode radiation with u < ui c . e 
propagates with a reduced refractive index. Adopting the notation 
j = (1, 2) to label mode (X,0) respectively, the transverse elliptic- 
ity may be expressed compactly: 



Ki = 



6+(-l) 8 -^agn(i/)v'6 a + C 



where 



2gr\ cos 9 
x — e 2 — g 2 — £r)(l — h/a) 



V = 



a (e sin 2 9 + rj cos 2 9) 
£ = 1 — h sin 2 9 j a 



(12) 

(13) 
(14) 
(15) 

(16) 



and the longitudinal polarization component has ellipticity given 
by: 



Ki = 



\ig + (e — rj)Ki cos #] sin 9 
e sin 2 9 + r\ cos 2 9 



(17) 



2.1 Cross sections 

The principal opacity sources in a fully ionized, light element 
plasma are the free-free absorption and ordinary Thomson scat- 
tering. The differential cross-section for scattering by electrons of 
photons with moment um k in p olarization mode i into momentum 
k! and mode j is JVentural 19791) : 
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(18) 



where ro = c 2 /(mc 2 ) is the classical electron radius (cm). Super- 
script (p) denotes plasma components only; e.g. the plasma polar- 
ization tensor is 



n 



(p) 



(19) 



The eigenvalues of <19> for the plasma component of the total di- 
electric are 



7T± 
7Tq 



= v- 1 [1 - 



e (P)±g(P)] 
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To excellent approximation, the squared eigenvalues can be writ- 
ten compactly as a rational expression in terms of the plasma con- 
stituents; including the damping term 7 S = v s /u) these are: 
2 

(20) 
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(E s m <=/ m »)^ 
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For coherent scattering, eqn i 1 81 is a sum over polarization com- 
ponents p, and has a scale fixed by the Thomson cross section 
ct t = 87tr 2 /3: 



do: 
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Integration of <2H over final momenta k! yields the partial cross- 
section to mode i photons which scatter into momentum k and 
mode j 



CTij(fc) = CTT 



hi 



p=-l 



ep(fe 
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(22) 



The bracketed term is the geometric correction to the scattering pro- 
cess induced by the magnetic field in the mode conserving (i = j) 
or mode crossing channels (i 7^ j). The damping factor 7 S is ob- 
tained from the sum of the radiative and collisional damping fre- 
quencies: 

v a , ra ,d = 2(ew) 2 /(3m s c 3 ) 



V s ,coU 



aov r ad/{n s (JT) 



where the free-free absorption coefficient ceo is defined below. 
Damping is only relevant near the resonant frequencies of tv 2 -, i.e., 

The total electron scattering cross-section to mode i photons 
is obtained by summing 1221 over final polarization j. 
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The magnetic free-free absorption coefficient is similarly derived; 
the total absorption for mode i photons is: 



£i i dn'£X>|4(*') 

j(v 3=1 p=-i 



2 I i /i \ I 2 2 
e p(^J ^P 



ooGf(fc) 



(24) 



where g is the Gaunt factor and ao is the absorption coefficient in 
the zero field limit, corrected for stimulated emission: 



q = n^riiZ^ 
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The corresponding emissivity is related to the absorption coeffi- 
cient through the Einstein rate coefficients; i.e., Kirchoff 's law: 

■i r,2 8e 6 / 2n X 1 / 2 -hv/kT n s n\ t~,c\ 

^ =n ^ Z S^{3^kf) e G ' {k) (26) 

Equation I26i is the correct thermal emissivity of unmagnetized 
plasma; for arbitrarily magnetized plasma, the expression must be 
multiplied by 0.5 to obtain the emissivity per mode. In the polar- 
izing medium, thermal emissivity is greater for the more opaque 
mode. It is convenient to scale the opacity and emissivity by the 
local mass density: 



off +n e <7i 
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The thermal emissivity j\ h is given by eqn 1261 . and the scattering 
emissivity in mode i is 



2 



dfi'N u'ik') 

3=1 



da 



fc'fc 



(29) 



where the integration is performed over incoming photon momenta. 

For unmagnetized plasma, the Gaunt factor s and their deri va- 
tives are evaluated from the fitting formulae of Itoh et al. 1 2000) or 
the tabulations of Sutherland 1 1998). The magnetic Gaunt factors 
and their temperature derivatives are evaluated by the integral ex- 
pressions found in Kirk & Meszaros ( 1980); this calculation pro- 
ceeds from an average of the ground state Coulomb matrix ele- 
ments over a one-dimensional electron distribution in the Landau 
ground level, which adequately describes the available electronic 
states in highly magnetized plasma. In weak magnetic fields, we 
expect that the anisotropy of the electron distribution is sufficiently 
mild in this regard to allow substitution of the Gaunt factors g± t o 
with the zero-field values. The transition from the classical to quan- 
tizing magnetic field regimes is characterized by the parameter 



I = kT eff /(huj c , e ) ~ 7.4 x WT eff /B 



(30) 



which is computed from the effective temperature T e g defined in 
j3|and neglecting variation in local temperature in a realistic at- 
mosphere. In this formulation the Gaunt factors are incompletely 
calculated over a range of (T e ff,B) corresponding to the transi- 
tion between the classical field and strongly magnetized regimes 

The geometric correction factors Gi(k) and Gf (k) in 1231 
and 03 

are fundamental quantities in the opacity calculations 
which fully describe all resonant behaviors of the magnetized 
plasma through self-consistent inclusion of ions and vacuum in 
the dielectric J2j, and which are computed explicitly for all model 
photon energies and trajectories. The cross sections are invariant 
to reversal of the magnetic field direction, and each cyclotron line 
has an effective width Alu C]S oc w c ,s in excess of either the ther- 
mal or natural width *y s of the resonance. Note that our choice for 
Il' p ) generalizes the polarization tensor for the multi-component 
plasma and produces the expected resonances at the electron and 
ion cyclotron frequencies. The vacuum dielectric influences only 
the calculation of the polarization vectors lilt with consequences 
discussed below. Cross sections for direct absorption and scatter- 
ing by ions are reduced from their respective electronic values by 
a factor 0(rn e /m p ) 2 and may be neglected without consequence 
when calculating the opacity. The opacity to radiation in paral- 
lel propagation to the magnetic field is approximately identical in 



each mode, regardless of the magnitude of the vacuum corrections 
to the plasma response; substantial differences in the differential 
cross sections can produce unequal intensities and finite polariza- 
tion along magnetic field lines. 

Electron and ion cyclotron absorption each result from differ- 
ent functional properties of the magneto-geometric factors. Both 
features are sensitive to projection of the X mode polarization vec- 
tors onto the circular basis about the magnetic field <1 It . In the 
electron feature, projection over all angles onto the resonant eigen- 
value -k 2 _ is complete; therefore only radiation in X mode is res- 
onantly absorbed while the O mode has no projection on irL and 
is consequently insensitive to resonant absorption. At frequencies 
uj <C ui c>e the §q is everywhere zero except near either ion cy- 
clotron frequency where it resonantly acquires a finite projection. 
Owing to the large relative opacity to waves polarized parallel to the 
magnetic field at these frequencies, it is this "polarization vector" 
resonance which dominates the total opacity to X mode radiation in 
the ion cyclotron resonance. The polarization structure in the res- 
onance is complicated further by the abrupt exchange of ± com - 
ponents at a critical frequency near ui C:P iBulik & Pavlovlll99d) . 
but the magnitude of the circular components remains symmetric 
across the feature. 

The magnetic field provides a preferred orientation to scatter- 
ing and absorption processes, resulting in anisotropic plasma re- 
sponse and finite polarization of propagating radiation. The mag- 
nitude of these effects scale principally with u e = {uj c ^/u]) 2 . Ion 
cyclotron features arise independent of vacuum corrections from a 
parametric resonance between electrons and ions in the dielectric 
tensor. Moreover, any additive contributions to the total plasma di- 
electric (e.g. line or ionization transitions) introduce additional po- 
larization resonances. The electronic and ionic plasma constituents 
cannot be consistently treated as independent entities. 

2.2 Vacuum resonance and mode ambiguity 

In the previous section, formation of the ion cyclotron resonance 
was shown to arise primarily from a critical parametric phe- 
nomenon in the spectrum of the polarization vector components, 
eP p . This is one example of the more generic competition between 
multiple sources of polarizability in magnetized plasma. More gen- 
erally, any additive component to the dielectric ^2} (e.g. transitions 
in neutral atoms) can produce other resonant behaviors, and intro- 
duce additional complexity to the spectrum of critical behavior. In 
intense magnetic fields (B > Bq), vacuum polarizability has a 
significant influence on the X-ray plasma opacities, and possesses 
parametric polarization resonances analogous to those of ion cy- 
clotron lines. Birefringence of the magne t ized v acuum has been 
investigated by many auth ors (e.g. lAdleJ ll97lh for review and 
IPavlov & Shibanovl {j.979) for discussion relevant to NS applica- 
tions); strictly, the effect remains finite for arbitrarily weak fields 
although, to excellent approximation, the vacuum corrections to 
weakly magnetized plasma are identical to the uncorrected (pure) 
plasma. Our discussion of the vacuum effect proceeds from the 
work of lSoffel et ail <1983h . 

Inspection of J12t reveals that the transverse ellipticities of 
the two modes are identical for arbitrary 8 when b = ±ivC — * 
Ki — ±&. This is distinct from (e.g.) the case of parallel propa- 
gation (k || B) in which the modes are circular, having identical 
magnitude, but opposite and unambiguous handedness. Where this 
occurs in a radiative transfer calculation, it may be impossible to 
uniquely specify the polarization content of the affected radiation 
field in two normal modes, sometimes called "mode collapse". To 
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understand the limitation of our model calculations, it is necessary 
to investigate the plasma conditions under which the mode defi- 
nition i 1 21 is rendered ambiguous. To deduce the critical density 
at which mode collapse occurs, it is sufficient to study the neces- 
sary condition that Re(fe) = for arbitrary propagation angles, 
i.e. where x — (eqns J13I14I ), To first order in the fine-structure 
constant, this condition is satisfied to good approximation when 



Vres ={q~h) (^^) 



(3D 



The critical density 13 1 i is undefined for lo > uj c , e and photon fre- 
quencies in this range are insensitive to vacuum resonance effects 
at all plasma densities. For photons ui < uj c ,e, the resonant density 
<3 1 1 delineates a transition from the plasma dominated dielectric 
response into the vacuum dominated regime. For polarization vec- 
tor calculations which invoke only real valued b, the modes acquire 
a finite but insubstantial measure of non-orthogonality when 6 = 0: 



= hsa?0/a < 1 



where the transverse components e 3 t = (T—kk)-e 3 . When comput- 
ing radiative transfer in such a medium, the mode identities remain 
distinct. Equation J3 li is a quadratic expression for two critical fre- 
quencies which are responsive to the vacuum resonance at a given 
plasma density: 



1±1 1- 



q — h \ ujc, 



(32) 



(q > h) for all B and the critical frequencies are defined only 
for LOp e < Lu1 e (q — h)/4. For v > v res the opacity is essen- 
tially unchanged from the pure plasma calculation and this regime 
is said to be plasma dominated. In the vacuum dominated regime 
(v < v re3 ) the modes become strictly linear, albeit slightly non- 
orthogonal, for all propagation angles 9 > 0, while both modes 
are fully circular for 9 = in all circumstances. An additional 
consequence of conversion to linear polarization is that the angular 
dependence to the opacity in each mode is quenched in the vacuum 
regime. This implies that radiation decoupled from the plasma in 
the regime v res > v(j v = 1) will also become isotropized but, as 
we find in i|5.3l this possibility is not realized in model atmosphere 
solutions. 

The imaginary part of b is approximately 



Im(fc) 



1 - 



sin 



2 cos 9 



(33) 



For magnetic field strengths capable of inducing the resonant phe- 
nomenon in NS atmospheres, the critical frequency lu+ ~ Lu c ,e 
is well above the limit of measurable thermal emission. Evaluat- 
ing 1331 at the low frequency resonance and taking Im(6) = ^JZ,, 
we find the critical angle at which the modes collapse at lj_ is 
cos 9- ~ 0.5 Im(6) <C 1, or 9- ~ 7t/2.Im(6) is a rapidly varying 
function near the critical angle, and for only a very narrow range of 
9 ~ will the polarization modes become indistinguishable. 

In the vacuum regime, the two critical frequencies are ui+ ~ 
Lu c ,e and w_ < w+. Both polarization modes become sensitive to 
the electron cyclotron resonance in the vacuum dominated regime; 
this is however of no special consequence for our NS spectra where 
io c ,e is unobservable for magnetic fields capable of producing the 
resonance behavior. From eqns 13 1 1 and I32i we find that, for a par- 
ticular magnetic field strength B, there are two critical frequencies, 



while, for any particular photon energy, the vacuum resonance core 
is found at a specific plasma density n e . 

Presently, no NS model atmosphere calculations fully specify 
the mixing of po larization at the mode collapse points (u)±,6±). 
iHo&Lail <2002) have pursued calculations invoking the two ex- 
treme limits of "no collapse" and of wholesale exchange of polar- 
ization identities through all photon trajectories, with the correct 
solution presumably intermediate to these extremes; these authors 
have since drawn conclusions similar t o those above w hich favor 
the limit of distinct polarization modes <Ho&Lail2 003). From the 
preceding discussion, it is acceptable to assume that mode-collapse 
will have at best a modest effect on the emergent NS spectrum re- 
gardless of the magnitude of the effect within the annulus, and 
we consider only the non-collapse limit in the remainder of this arti- 
cle. The mode ambiguity is related to the breakdown of the Faraday 
depolarizing limit, and a more exacting solution will likely require 
radiative transport of the full Stokes vector, and careful treatment 
of the anti-hermitian components of the plasma dielectric. 



3 RADIATIVE TRANSFER AND NEUTRON STAR 
ATMOSPHERES 

The fundamental parameters of the model atmosphere are the sur- 
face gravity g 3 , the total flux propagating through the atmosphere 
J- eons, the magnetic field strength and orientation B, and the 
plasma composition. For a given equation of state, these parame- 
ters uniquely specify the model. 

The surface gravity of a relativistic star is 



: 7- 



GM 



GM/R 2 



R2 y/1 - 2GM/R 



and for M, R characteristic of NSs, g s ~ O(10 14 ) cm s~ 2 . The 
large surface gravity gives an atmospheric scale height d <C R, and 
gradients in the gravitational and magnetic fields of the star may be 
neglected for the purpose of calculating the atmospheric structure. 
The thin plasma atmosphere is well described by a semi-infinite 
plane parallel geometry. We do not consider extended atmospheres 
or those with bulk motions. 

In the absence of heat sources or sinks, the total energy flux is 
a conserved property, and we define the effective temperature T e g 
in terms of this conserved flux: 

PSBTejj = J- cons 

for a sb the Stefan-Boltzmann constant. Conductive transport at 
densities characteristic of NS atmospheres (p ~ 0.1 — 10 g cm" 3 ) 
is inefficient compared to radiative diffusion, while convective in- 
stabilities, already suppressed in the limit of complete ionization, 
are further inhibited by Ohmic dissipation in strong magnetic fields. 
Consequently, neither conductive nor convective transport is in- 
cluded in our calculations. 

In a semi-infinite, plane parallel geometry, the steady-state 
transport of radiation in mode j for photon momentum k at any 
particular depth in the stellar medium is described by the differen- 
tial equation 



(34) 



The momentum vector k is assumed to be outward pointing, and is 
specified by the surface latitude p = cos 9k and azimuthal angle cj>; 
the latter is measured with respect to a reference direction to be de- 
fined later. All functional dependence on the magnetic field vector 
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is subsumed in the formulation of xi(k) and rji(k). It is advanta- 
geous to write the radiative transfer equation 1341 under transfor- 
mation to the pair of conjugate intensity (Feautrier) variables: 



nt(k) = ~ [ll(k) + Il(-k)] 

vi(k) = l[ii( k )-iu-k)] 

where for k = k(fj,,(f>), —k — k(—fi,n 
transfer equations in these variables are 

dvi(k) 



/' 



/'- 



dm 
dul(k) 
dm 



= X l(k)ul(k)- V l(k) 
= xi(k)vi(k) 



(35) 
(36) 

and /j, > 0. The 

(37) 
(38) 



where we have relied on the residual symmetry of the plasma emis- 
sivity about the magnetic field orientation, i.e. r]l(k) — r)i( — k). 
Following the usual Featrier development, we obtain the transfer 
equations for the surface boundary, intermediate depth, and lower 
boundary, respectively: 



-l duj _ j 
dm ~ U » 



dm 



dul 



-i duj 

X] ' v dm 

3 cons 



1 1 



-XR- 



dB v 



dm 16 T 3 ~" dT 
for the Rosseland mean opacity: 




XR 4(tsbT 3 



(39) 



(40) 



(41) 



(42) 



At the stellar surface, the term ll( — k) can be used to quan- 
tify either the intensity in mode j of the incident radiation field of 
an externally illuminated atmosphere, or a small correction to the 
intensity field to account for the finite optical depth of the surface 
stratum in a numerical calculation (0 < rl <C 1). For our present 
discussion, we will set Il(—k) = 0. 

The large optical depth in either mode to the boundary allows 
for ample redistribution of the flux by (principally) thermal absorp- 
tion and reemission; in magnetized models, the boundary flux of 
<41> can be arbitrarily distributed in the two modes with no ef- 
fect on the result of the radiative transfer. In strongly magnetized 
plasma, the opacity introduces an additional non-trivial angular de- 
pendence to xr m 1411 . neglected at the lower boundary but other- 
wise accounted for in the RTEs at all other depths. The anisotropy 
ultimately has no appreciable effect on the emergent radiation field, 
owing to the tremendous optical depth to the lower boundary (typ- 
ically r > 10 2 ), and the plasma can be modeled as isotropic with 
no effect to the emergent spectrum. 

3.1 Radiative equilibrium 

The total flux gradient is found by integration of the transfer equa- 
tion J34I over all frequencies and angles where the opacity and 
emissivity have been expanded into their respective thermal and 
scattering contributions. Using the definition of the (coherent) scat- 
tering emissivity jf c v , and demanding that the flux be invariant with 
depth, this is 



dfi 



2 



■ th 



= 



(43) 



Equation 143 > is a necessary condition for radiative equilibrium, 
but is not by itself sufficient to constrain the magnitude of the flux 
— the thermal balance prescribed here does not relate continuous 
variations in flux with depth. When simultaneously integrated with 
the system of radiative transfer equations (RTEs), where the total 
conserved flux is specified in the depth boundary 14 II . equation 
J43I provides a robust constraint. To directly calculate the flux T v 
we evaluate the conjugate field vj„ from either the Feautrier defi- 
nitions at the surface or from eqn J38I at any other depth: 



/>27l /*1 2 

■Tv = 2 / / vyy^v dju 
Jo Jo =1 



(44) 



3.2 Hydrostatic equilibrium 

Our model of the neutron star atmosphere will neglect bulk motion 
of the plasma, and will therefore assume a hydrostatic structure. 
The total pressure at any depth of the atmosphere is the sum of con- 
tributions from the ideal gas pressure, the radiation field, and non- 
ideal effects owing to interactions between plasma constituents: 



P — Pideal ~\~ Priori 



~t~ Pra 



(45) 



The ideal gas pressure, including electron degeneracy, is 

Pideal = NkT + n e (y p - l)kT (46) 

where N is the total number density of the plasma constituents, 
and the parameter y p is the ratio of Fermi integrals I P (a), y P = 
I p / (pl p -i). For unmagnetized plasma, p = 3/2, and in quantizing 
magnetic fields, p = 1/2. The electron density is an independent 
model variable in our calculation; to ensure a self-consistent equa- 
tion of state, we compare th e value of n e of the current iterate to 
that given by Jciavto jl98llPotekhin et alll 999) 



4 T 



1/2 



l> 1 
i < 1 



(47) 



where A e = (2nf3h 2 /mekT) 1 ^ 2 is the electron thermal wave- 
length, and the parameter £ is defined in <30t . The electron de- 
generacy parameter a is evaluated from n e either by inversion 
of the Fermi integral l Antia 1993) or by a numerical root search, 
from which the degeneracy pressure is calculated. In strongly quan- 
tizing magnetic fields, the phase space volume occupied by the 
electron distribution is much smaller than that of a weakly mag- 
netized plasma, tempering the onset of degeneracy which occurs 
at higher densities than for unmagnetized plasma. At intermedi- 
ate magnetic field strengths, the degeneracy contribution to the gas 
pressure should be evaluated by a sum over occupied Landau states. 
In practice, degeneracy pressure contributes less than ~ 4% at even 
the deepest model strata, and has no discernible influence on the 
emergent spectrum. In our models the degeneracy pressure at finite 
B < 10 10 G is approximated by the zero field limit. 

Non-ideal contributions to the pressure include the Coulomb 
interaction of the ionized plasma and, for an incompletely ion- 
ized plasma, the atom-atom and atom-ion interactions in the neutral 
plasma component. The latter are neglected in our simplified EOS; 
for the former, we use the Debye-Hiickel approximation in both 
magnetized and unmagnetized plasma lClavtonll983l) : 

2e 3 / 2n\ 1 / 2 3/2 



(48) 



where £ is a compositional parameter 0(1). This approximation 
breaks down when the Coulomb energy at the mean electron-ion 



Model atmospheres and thermal spectra of magnetized neutron stars 1 



separation becomes comparable to the local thermal energy, which 
can occur at depth in models which are both strongly magnetized 
and very cool. 

The radiation pressure is formally a tensor quantity but usu- 
ally taken to be isotropic in conventional atmosphere models. Fi- 
nite anisotropies to the pressure tensor are induced in oblique mag- 
netic fields; however, as we shall see from the discussion of model 
results, the radiation force plays only a minor role in the total pres- 
sure structure of the neutron star atmosphere (this is especially true 
in the region of photon decoupling) and we will consider only the 
zz component of the pressure tensor regardless of field orientation 
in these calculations. In CGS units this scalar pressure is: 



Prad 



An 
c 



dv 



2 



+ S2 



dQ jj 2 \ ujj 



(49) 



The fiducial depth in our atmosphere model is the Lagrange 
variable m (g cm - ), in terms of which the total pressure gradient 
is: 



dP 

dm 



(50) 



Eqn I50i is to be integrated simultaneously with the RTE system 
and other constraints. In the limit of m — > the total pressure will 
equal that of the emergent radiation field; in numerical calculations 
it is not possible to realize this limit, and the pressure P(mo) will 
retain finite contributions from each source, dominated by the gas 
pressure P g . At the surface, we require P(mo) = <? s mo. 

Finally, the system of equations is closed by the EOS. In the 
limit of complete ionization, the statistical equilibrium for a light 
element plasma is particularly simple: we require only that the local 
charge distribution of the model be neutral, and that particle number 
be conserved everywhere: 



N 



Although convective transport is neglected in our prescription 
for energy transport and flux conservation in the stellar atmosphere, 
conditions amenable to convection may arise in the resultant atmo- 
spheric structure, and the validity of our non-convective model hy- 
pothesis should be evaluated in the final result. The onset of convec- 
tive instability is estimated by the Schwarzschild criterion, gener- 
alized to include the effects of radiation pressure, and of the neutral 
species gradient in unmagnetized atmospheres: 



VR ^{d^p) R > \d^p) A ^ VA 



(51) 



The radiative gradient Vr is evaluated from the model T(P) 
by numeric differentiation and Va is calculated as described by 
iKrishna-Swamvl ll96ll) . This measure of the adiabatic gradient 
does not include the effect of Ohmic dissipation in magnetized 
plasma, which tends to improve the stability of the plasma against 
the onset of convection. None of th e models pre sented in j5|show 
evidence of convective instability; IZavlin et alJ ll996ah find that 
convection occurs in light element atmospheres only for exception- 
ally cool models, T eff < 5 x 10 4 K. 



4 SOLUTION TECHNIQUE 

Apart from the strategy used to refine the model, successful conver- 
gence will often depend critically upon the properties of the trial or 



seed solution, and on the distribution of points on the discrete grids 
over continuous model dimensions, e.g. energy, angle and depth. 

4.1 Seed thermal structure 

We use a power-law prescription for the plasma conductivity as- 
suming free-free absorption is the dominant source of opacity, i.e. 
the Kramers opacity l Clavton 1983). This simplified model has two 
functions: first, to suggest reasonable bounds for the depth variable 
m which will cover the range of photospheric depth for a prescribed 
range of photon energies; and second, to construct an initial tem- 
perature profile which crudely approximates radiative equilibrium. 
Derivations of the mass, temperature and optical depth relations 
for this model proceed from the discussion in iHevl & Hemauistl 

For the monochromatic optical depth r to photons of energy E 
(keV), the equivalent mass column density (g cm~ ) in this model 
is given by 



m E (r) = 2.74 x ICT 2 ^ 7 (E 3 t) 17/28 



(Ms 



-1/2 



(52) 



Here, g s ,i4 = g s /W 14 and /i is the mean molecular weight of 
the plasma, which can be defined in some average sense for H/He 
mixtures with compositional gradients. The required depth range is 
therefore approximately bounded by the range of photon energies 
considered in the model. For the NS spectra under consideration, 
we take for this energy range 1CP 3 — 10 keV. The Kramers model 
underestimates the surface temperature of the converged result, and 
generally overestimates the opacity of the plasma in the outermost 
strata. Consequently, it is sufficient to define the surface boundary 
by evaluating 1521 for the lowest model energy at r = 0.01, and 
the lower boundary by evaluating 1521 for the highest model en- 
ergy at r = 1.0. For models with log_B 11.0 the lower depth 
boundary is multiplied by Bn = 6/10 G to account for the re- 
duced opacity at X-ray energies in the extraordinary polarization 
mode. The total optical depth range for a converged model is typ- 
ically t ~ 10~ 5 — 10 5 . The temperature profile derived from the 
power-law model is 

T(m)~ 56 (fig s T? ff m 2 ) 2/17 (53) 

The initial density profiles can be approximated by assuming 
only ideal gas pressure contributes to the total. 

4.2 Mesh generation 

The rate and quality of model convergence for a particular choice of 
parameters are sensitive to the choice of discretization of the con- 
tinuous model dimensions, and will often require the careful distri- 
bution of grid points in depth rn, photon energies E and trajecto- 
ries through the plasma (/i, (f>). The current numerical implemen- 
tation requires the angle-energy mesh to remain fixed during the 
calculation. We expect that a reliable solution should result from a 
mesh which is designed to capture significant features of the plasma 
opacity with respect to these dimensions. 



4.2.1 Frequency mesh 

For most NS spectra we consider the energy range E = 10~ 3 - 10 
keV which provides sufficient coverage of the spectrum from op- 
tical through soft X-ray energies for a broad range of T e g. The 
continuum frequency mesh is allocated by dividing the total en- 
ergy range into logarithmic subintervals, each populated with equal 
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numbers of mesh points. Any of these intervals may be further sub- 
divided to accommodate an atomic or cyclotron feature, in which 
case the continuum interval is split into smaller intervals with mesh 
point densities proportional to their widths. We typically use 6 — 8 
points per logE to describe the continuum, with an additional 
NC = 4 — 8 points cast in each cyclotron resonance. 

4.2.2 Angle mesh 

Several strategies for prescribing the angle mesh may be consid- 
ered, and we take advantage of symmetries about the magnetic field 
when available. The radiative transfer equations in terms of the con- 
jugate Feautrier variables require definition of the outgoing photon 
trajectories p = (0 — 1]. 

For weakly magnetized plasma (log B < 10), the anisotropy 
imposed on the emergent radiation field is mild at X-ray energies. 
In this nearly isotropic plasma, we adopt a simple Gauss-Legendre 
or Gauss-Radau distribution in the viewing latitude p,; the latter 
distribution is constructed by holding the a = 1 node fixed in 
the Legendre polynomial. Model results calculated in the isotropic 
limit which include the normal LOS do not differ substantially from 
those where this trajectory is excluded. In this regime, as few as 5 
points can be used to describe the angular intensity distribution. 

Departures from the isotropic plasma response affect the X- 
ray spectra when log B > 10, where large variations in the angular 
dependence of the opacity allow for photon decoupling from broad 
range of plasma temperature. For u e > 5.3 the geometric factor 
Gx acquires a maximum at 9 P roughly delineating a bimodal an- 
gular distribution of radiation: a pencil beam component directed 
along the magnetic field lines, and a broader fan component at more 
oblique angles. The angular width of the pencil component declines 
with increasing u e after achieving a maximum near u e ~ 11; res- 
onances in the polarization vectors (ionic or vacuum) also have 
large opening angles, but the pencil beam transports little flux in 
NS models with B > 10 13 G. It is desirable to prescribe a fixed 
angle distribution which samples the dependence across the entire 
spectrum of model energies for a given magnetic field strength, 
and which emphasizes the transition between the two components. 
These pencil weighted schemes distribute mesh points on intervals 
of log (1 — p) where (1 — p) m in is generally 10~ 4 . Good reso- 
lution of the pencil and fan emission components can be achieved 
with as few as 10-16 points. As in the isotropic case, the normal 
sighting can be included as desired. 

Oblique fields & the azimuthal mesh — For normal magnetic 
field vectors, the radiation field is azimuthally symmetric models 
and degenerate in the cf> dimension; for this, we use a single mesh 
point at cos cj> = 1. Weakly magnetized models can also be com- 
puted with sufficient accuracy in this approximation. Oblique mag- 
netic fields B > 10 10 G clearly break the azimuthal symmetry and 
for generally oblique magnetic field inclinations we are obliged to 
consider <j> : — 7t/2 with emphasis on the pencil width about 
the field axis (8t, <t> = 0). When 8t — n/2, the plasma response 
recovers partial symmetry, and the calculation can be restricted to 
4> : — n/4, with the remainder of the plane mapping into this 
range. The direction <j> — is defined by the projection of the 
(oblique) magnetic field vector onto the plane of the atmosphere 
(x — b — n cos 0b). For generally oblique magnetic fields, we use 
5 — 8 azimuth points in the following overlapping distributions: (i) 
two nodes are assigned at cos</!> = —1,1 and (ii) the remainder 
distributed logarithmically in the quadrant cos</> = [0, 1] similarly 
to the pencil weighted scheme described above. Fully orthogonal 
models (b — n/2) can be calculated with as few as 4 azimuthal 



points, and for these we choose a Gauss-Lobotto distribution with 
anchored nodes at cos <f> = [0,1]. Regardless of their distribution, 
we calculate the angle the trajectory and the magnetic field for each 
angle pair (p, <p) : 



cos 9bk = cos 4> (sin 6^ sin 6k + cos 8^ 



(54) 



for use when evaluating the geometric corrections to the plasma 
opacity detailed in J3| Model results are invariant under reversal of 
the field direction. 



4.2.3 Depth mesh 

In general, we compute the atmosphere models on a fixed depth 
mesh with 10 points per decade for a total of 90 — 130 points. 
Models including a strong magnetic field, with non-trivial vacuum 
contributions to the opacity, often require special care; experience 
has shown that these models benefit from an additional depth mesh 
component which adaptively traces the coupling of density with 
photo n energy n ear the vacuum resonance (e.g. eqn 13 1 i . see also 
iHo & Lail 12003)), These adaptive depth points are found by inter- 
polating the current density profile n e to find the critical density for 
vacuum resonance formation at each model energy from fm . and 
ignoring photon energies which do not see the resonance on the cur- 
rent range of n e described by the model. A single adaptive point is 
assigned to the "core" of the vacuum resonance; additional points 
can be symmetrically distributed at p res (l ± ne), n = 1, 2 . . ., 
for e <C 1 but this is not required. All independent model vari- 
ables (T, n, u) are then interpolated for each adaptive depth point 
from the current model solution. In successive iterates, the previ- 
ous adaptive mesh is discarded after interpolation of the current 
iterate's adaptive points. 



4.3 Formal radiative transfer calculation 

For a given atmospheric structure T(m),n a (m) we calculate the 
opacities and thermal emissivities at each strata of the model. The 
radiative transfer equations 13914 It are transcribed to a simple dif- 
ferencing scheme; the discretized RTE at intermediate depth for 
radiation of a particular energy, mode and trajectory (with these 
indices suppressed) is: 



-p. D-Ud-i + \p {D- + D+) + Amx] u d — p D + u d+1 
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(55) 



and the surface and depth boundary equations are 
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(56) 



(57) 



dT ' 2" 

for d — 2 . . . ND — 1. The differenced boundary equations J561 
I57i have been promoted to second order accuracy by substitution of 
<40> ; higher order accu racy can be obtained by evaluating the gra- 
dient of the emissivity lAueJl97d) . but does not substantially im- 
prove our results. The depth intervals are Am± — \md±i — 
Am = Am+ + Am_ and the differencing coefficients D are 
lAuer et all 19771) : 
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D± 



Am± 

The scattering integrals are expanded as a sum over incoming 
photon modes and trajectories which explicitly couples radiation in 
the two polarizations states, and integrates the scattering emissivity 
in-situ: 



p v 



2 

£■ 

3=1 



cK2' 



The angular integration weights are the product of the respec- 
tive /j, and <f> Gaussian weights corresponding to the trajectories k. 
Rendering integrations of the total opacity and emissivities on the 
fixed angle mesh accurately preserves both thermal balance and the 
conservative scattering processes. 

For a particular model energy E, the coefficients of the vec- 



tor Ud,E = [W 
coefficients 



: e from J55I57I are written into matrix 



Ad,EUd-l,E + Bd,EUd,E — Cd,EUd+l,E — Kd,l 



(58) 



which are generally full but diagonally dominant; the Kd are the 
thermal emissivity factors and boundary terms from eqns 1561 and 
1571 . For coherent scattering, radiation streams of different ener- 
gies are not coupled and, per frequency, the RTE system is simply 
expressed in this block tridiagonal form for mode-trajectory pairs 
organized by depth which can be solved by standard elimination 
techniques. Matrix inverses are never explicitly evaluated - instead 
the common matrix terms in the projection operators are LU fac- 
tored and the intermediate vectors and matrix operators then found 
by back-substitution. All algebraic reductions in our code are made 
using standard BLAS and LAPACK routines. 

The resultant radiation field, while consistent with the pre- 
scribed atmospheric structure, cannot be expected to satisfy flux 
conservation. The trial solution will therefore be iteratively refined 
by deducing corrections to each model variable such that the entire 
system of RTEs and constraints are driven towards a convergent 
solution. This proceeds by complete linearization (CL). 

4.4 CL algebra 

The tridiagonal organization of the formal RTE system in the 
previous section provides the algebraic basis for extension to 
the full iterative problem in the complete-linearization technique 
lAuer & Mihalasiri969l) . although the detailed organization of the 
coefficients is distinct from that of the formal transfer calculation. 
The constrained RTE system is again transcribed to a numeric grid 
in a simple differencing scheme, but now linear perturbations to 
the radiation components and physical fields are the independent 
model variables. The correction to any model variable will be asso- 
ciated with a linearized equation: either an RTE (for u) or con- 
straint (for T,N,n s ). The solution is represented by the vector 
* = (*i . . . *£>), where * d = (JV, T, n Sj w^„)d are formed from 
the current model state. Corrections to the model solution are found 
by solving an equation of the form OS^f — E where E is the error 
in the current solution, and O is an operator derived from coeffi- 
cients of the perturbed equation s. For brevity we sketch the main 
elements of the technique; see Mihala^ J 19781) for a complete ex- 
ample of complete linearization applied to a simple atmosphere. 

In each model equation, the independent variables are re- 
placed by first order corrections, T — > T + ST, etc. Functions 
are linearized by taking partial derivatives with respect to the inde- 
pendent model variables: 



dn s 



(59) 



and similarly for the emissivities rf and rj sc . The efficiency of the 
CL cycle is diminished by properties held constant during the it- 
eration, therefore differentiation like that in J59I includes partials 
taken with respect to the squared polarization components and the 
Gaunt factors. All derivatives in the model are evaluated analyt- 
ically at the current values of the model variables; for numerical 
work these are most conveniently expressed in logarithmic form. 
Upon substitution of all model variables, the linearized equations 
are reorganized by collecting terms from the unperturbed equa- 
tions to form the error E, while coefficients of each model vari- 
able correction are grouped by depth. The CL coefficients of the 
constrained system of equations are written into matrix coefficients 



Ad5t>d-i + B d 8^d - CdS^d+i = E d 



(60) 



where the matrix coefficients A,B,C now contain the full spec- 
trum of differencing coefficients, integration weights and constants 
of the linearized RTE system and its constraints per depth, and Ed 
is the error vector whose components are given by the difference 
between the RTEs (constraints) for u (N, T, n s ) and the exact so- 
lution calculated from the current model state. Coefficients of the 
jth variable at depth d — 1 in the ith equation at depth d are written 
into the Ad,ij etc. The entire system forms a block tridiagonal ma- 
trix which is solved by LU factorization and substitution. The new 
model state is found by updating each model variable with its lin- 
ear corrections, and the system is iteratively corrected in this fash- 
ion until the root vector is rendered within some prescribed toler- 
ance to be described below. This procedure is multivariate Newton- 
Raphson iteration. Note that information (model corrections) from 
a given stratum propagate over the mean free path of the radiation 
field via the transfer equations, and model convergence is global. 

CL often over-steps the convergent solution in early stages of 
the computation, and we limit the magnitude of the temperature 
corrections to \AT/T\ ^ (0.4 — 0.7) to prevent model variables 
from acquiring non-physical values from "rogue" iterates. When 
this condition is engaged at any depth, corrections to all other 
fields are discarded, the density made consistent with g s and the 
correction-limited T(m) as in £14,11 followed by a formal recalcu- 
lation of the radiation field ( i)4.3> . The typical CL cycle is sum- 
marized by (1) allocation of mesh points in all model dimensions, 
either from command-line arguments to the program, or through 
defaults; (2) this is followed by tabulation of the Gaunt factors and 
their temperature derivatives for each model photon energy in mag- 
netized models, over a range of temperature log T e g ± 1.5; (3) cal- 
culation of the trial thermal structure T(m),n s (m), and the run 
of opacity xL an d emissivity rjl with depth for all photon trajec- 
tories 6bk in each mode; (4) a formal calculation of the radiation 
field from these opacities ( £|4.3> ; followed by (5) the CL iterates. 
We adopt as acceptable convergence criteria that AT/T < 10 -3 
and AF/F cons < 10 -3 at all depths. These limits are achieved in 
5-13 iterates for typical choices of computational grid densities and 
dimensions, while the convergence rate can become superlinear in 
the final iterates. 



5 MODEL RESULTS 

We have implemented the methodology of [j4]in an efficient stellar 
atmosphere code, and computed model atmospheres and spectra 
for a range of parameters consistent with measured and inferred 
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properties of cooling, isolated NS. We consider stars with surface 
magnetic fields up to 10 14 G, and assume an ionized Hydrogen 
plasma photosphere described by the opacities formulated in j2l 
For weak and itermediate field strengths, we consider the range 
5.5 ^ log Teff ^ 6.5, while for the magnetar case we consider 
only 6.2 ^ log T B g ^ 6.7. We restrict our discussion to mod- 
els with surface gravity logg 3 = 14.38, corresponding to NS of 
canonical mass M = 1.4 Mq and radius 7? = 10 km. Spectra 
are illustrated in the stellar surface frame and must be redshifted 
by 2 = 0.306 to describe the flux measured by a distant observer. 
All results are derived from models converged to the tolerance de- 
scribed in i|4.4l No smoothing is applied to figures derived from 
converged model results. 

The principal energy scaling of the geometric corrections 
Gi(k) and Gf (k) enter through the magnetic field strength in the 
ratio u e = {uJc.e/u)) 2 , and it is natural to organize our discussion 
of model results by ranges of B. This grouping emphasizes the 
discriminatory properties of the stellar spectra, but also naturally 
collects results appropriate for comparison to different classes of 
observed sources (quiescent X-ray binaries, radio pulsars, anoma- 
lous X-ray pulsars). The selection of models presented in this dis- 
cussion, while not exhaustive, is representative of typical values for 
observed manifestations of NS thermal X-ray emission. 

X-ray polarimetry is a potentially predictive measure of the 
surface properties of cooling NS, but is of immediate interest only 
for relating the model input physics to the theoretical spectra. In 
our present study, it is adequate to define the net polarization II(iJ) 
in terms of the partial flux in each polarization mode I44i : 

1 ; T X {E)+T°{E) 

The flux input to the model J4H can have arbitrary polarization con- 
tent; the radiation field is efficiently redistributed by absorption and 
scattering processes and the final spectrum is independent of the 
input polarization fraction provided r^ asc 3> 1. The polarization 
content of the intensities at large depth is essentially zero, and that 
of the final emergent spectrum is determined by the relative tem- 
peratures at which the two modes decouple from the atmosphere. 
Upon decoupling from the plasma, the intensity and its angular dis- 
tribution (i.e. flux and pressure) in either mode asymptotically ac- 
quire values which are unchanging in the higher strata. As the gas 
pressure declines toward the atmospheric boundary, the radiation 
pressure assumes the dominant role; exterior to the star, the pres- 
sure is entirely that of the radiation field. All model calculations 
described here account for the radiation pressure, but throughout 
the range of photon decoupling its contribution to the total P is 
negligible. 

5.1 Unmagnetized models (B < 10 10 G) 

The onset of magnetic corrections to the plasma opacity occurs 
where u e — 1, and the electron cyclotron feature clearly delin- 
eates the polarizing properties of the magnetized plasma. Conse- 
quently, the X-ray spectra of magnetized NS atmospheres having 
B < 10 10 G are essentially indistinguishable from the zero-field 
case for the same plasma composition and total flux, and may be 
classified as unmagnetized models. The transition from classical to 
quantizing magnetic field occurs within the range of T e ff described 
for B = 10 9 — 10 10 G. As described in S0our calculation of the 
bremsstrahlung Gaunt factor is somewhat inaccurate for I ~ 1, and 
here we will assume the zero-field g for each effective temperature 
at these field strengths. Consequently, the polarization properties 



derived for models in this regime are influenced only by the nor- 
mal mode vector decomposition, but permit comparison between 
models at different Teff having the same field strength and guar- 
antee consistent calculation of g at all strata of a given model. In 
Figure Q we compare spectra calculated in the weak field regime 
for a range of T^and find no magnetic field dependence for the 
emission E > 0.2 keV. For the range of T e g considered here, the 
spectra in the unmagnetized field regime satisfy a linear relation be- 
tween the peak spectral energy the effective temperature similar to 
the Wien displacement of the Planck function; E p <, a k — i.7kT e ff 
in these unmagnetized models. 

Departures from the zero field results at optical and UV wave- 
lengths for finite B arise from the redistribution of resonantly ab- 
sorbed radiation to longer wavelengths in accordance with flux con- 
stancy in radiative equilibrium. The magnitude of this redistribution 
is greater where the cyclotron line is approximately coincident with 
the peak thermal emission (e.g. log T e g — 5.6 at 10 G) owing 
to both the scaling of the effective cyclotron line width and be- 
cause there is more radiation in the spectral peak which can be 
reprocessed. Additional details of the influence of cyclotron en- 
ergy redistribution on visual magnitudes are presented elsewhere 
(Llovd et al. 2003). More generally, photons absorbed in a neutral 
plasma component is likewise reemitted at longer wavelengths in 
radiative equilibrium producing enhanced emission below the ab- 
sorption feature. 

Only X mode radiation is resonantly absorbed in the elec- 
tron cyclotron line, and the excess ordinary mode contribution 
yields a large net polarization in the resonance. The effective width 
of the cyclotron feature is approximately scale invariant. At op- 
tical and near UV wavelengths the X mode photosphere is dis- 
placed to higher temperatures, owing to the reduced opacity in this 
mode. The spectrum is largely unpolarized at energies lo U) c ,e, 
for which the opacity reverts to the anticipated unmagnetized re- 
sult, although, as the magnetic field strength approaches 10 G 
(hoj Cl e ~ 0.12 keV), the unpolarized limit is recovered only grad- 
ually at high energies. For energies E > 2.5 keV, the photosphere 
temperature is saturated at T ~ 2 x 10 6 K where (frequency in- 
dependent) electron scattering is the principal opacity source. The 
polarized spectra and plasma decoupling properties of the 10 10 G 
models are illustrated for several T^in Figs.|2|and|3| The decou- 
pling density p(te) = 1 in each mode is approximately indepen- 
dent of flux temperature for a given value of B. Collective phenom- 
ena become significant in the plasma response at long wavelengths, 
are not fully described by the opacity formulated of [0- this limita- 
tion can be seen at optical frequencies in the 10 10 G models where 
photons decouple at v e ~ 1 (Fig.[3j- 

Radiative equilibrium mitigates the influence of cyclotron ab- 
sorption by adjusting the asymptotic surface structure T(rn). Fig- 
ure [4] demonstrates the magnitude of this effect for two values of 
T e ff. The ordering of surface temperature T(mo) with magnetic 
field differs for each value of T e g , and the departure from the un- 
magnetized temperature profile at finite B occurs at greater depth 
with increasing magnetic field. At depth where only photons with 
lu > u e , c decouple, the models have identical thermal structure. 
Radiation emerging from this isotropic plasma has an angular dis- 
tribution characteristic of limb-darkening, which is preserved for 
oblique magnetic fields of modest strenght. 



5.2 Intermediate fields (10 11 < B < 10 13 G) 

In the unmagnetized plasma regime, geometric corrections to the 
opacity are both weak in amplitude, and nearly independent of 



Model atmospheres and thermal spectra of magnetized neutron stars 1 1 



propagation angle; consequently, high energy emission in weakly 
magnetized models is well described by the equivalent zero field 
results. Departures from isotropic plasma response are prominent 
in the X-ray spectra of models with 10 11 < B < 10 13 G com- 
parable e.g. to the radio pulsars (RPs) where they are manifest in 
the finite polarization H(E) of the spectrum and in the angular dis- 
tribution of intensities. Polarization of the thermal radiation is es- 
sentially complete for B > 10 12 G, and the spectral envelopes 
of these models differ substantially from their weakly magnetized 
counterparts, owing to the harder energy dependence of the X mode 
opacity (oc v~ v ) which suppresses the amplitude of emission in the 
high energy tail and shifts the peak emission energy to lower fre- 
quencies. Inflections in the equilibrium thermal structure (Fig.0 
are related to displacement of X mode decoupling to deeper, hotter 
strata than for the ordinary mode. In this regime, the asymptotic 
surface temperature declines for increasing B owing to the reduced 
energy density in the O mode field. The peak energy displacement 
with effective temperature is E pea k ~ 4.0fc T e g , and the spectrum 
is dominated by radiation component is the X mode field. Vacuum 
corrections to the plasma dielectric are negligible in these models. 

The proton cyclotron resonance enters the optical regime at 
logB ~ 11.5, and soft X-rays at B ~ 10 13 G. The apparant 
depth of the feature is greater than that of the electron feature seen 
in weakly magnetized models owing to the overall depletion of O 
mode intensity at u < Lu c .,e- The proton line width scales with the 
cyclotron frequency similar to that of the electron resonance (Fig. 
|6j. Collective phenomena play an increasingly important role at 
long wavelengths in these models, and the energy below which the 
plasma attenuation becomes significant is E p i as — 0.5 log B — 7.6. 
At these long wavelengths the specific photosphere is coincident 
with the boundary v e owing to the large optical depth accumulated 
at this boundary. The broad electron resonance E ~ 11.6 keV for 
models of 10 12 G is not resolved in these calculations, though its 
influence is obvious in the polarized spectra in Figures[5HH 

The angular intensity distributions in the intermediate field 
regime acquire a pronounced pencil beam component oriented with 
the magnetic field vector. The pencil is essentially unpolarized due 
to the degenerate opacity for dt ~ but acquires large net polar- 
ization away from this trajectory. The amplitude of the pencil com- 
ponent declines with increasing magnetic field inclination due to 
previously mentioned limb effects, but the pencil width is approx- 
imately independent of inclination (Fig. |8}. The O mode opacity 
is not magnetically depressed outside the pencil and the intensity 
profile acquires a highly polarized fan-like profile dominated by X 
mode emission (Fig. [5}- Despite the large intensity in the pencil, 
less than 10% of the total flux is transported in this component at 
10 12 G; and by 10 13 G, less than 1% is emitted in the pencil. In 
general, the O mode flux contribution declines with increasing field 
inclination b, as does the total emission in the high energy tail (Fig. 

5.3 Intense field limit (B ~ 10 14 G) 

The abrupt X mode opacity variation induced by vacuum polar- 
ization plays a crucial role in the equilibrium thermal structure 
and spectral properties of ultramagnetized plasma atmospheres 
(B > 10 14 G). The magnetic suppression of the X mode opacities 
is greater in this regime than for the RPs, and O mode intensities 
contribute < 1% of the total flux. Strictly, the precise polariza- 
tion content of the radiation field is indeterminate in the normal- 
mode analysis as radiation propagates from the plasma dominated 
regime, through the resonance into the vacuum dominated regime. 



The shape of the spectral envelope is substantially affected by both 
direct absorption within the opaque layer, and its indirect influence 
on the equilibrium thermal structure. The displacement rules de- 
scribed for B < 10 13 G is ambiguous for 10 14 G owing to spectral 
reshaping through vacuum effects. The pencil component is of neg- 
ligible size in magnetar X-ray spectra and contributes <C 1% of the 
total flux. 

Energy absorbed from the X mode field in the vacuum reso- 
nance heats the local plasma. There is substantial variation in the 
precise depth at which the absorption will occur across the soft 
X-ray band 13 1 1 and the net effect is the formation of an opaque 
"heating layer" in the atmosphere (Figure fTTl . The influence of this 
absorption on the emergent spectrum depends on the proximity of 
the absorbing layer to the X mode decoupling layer. Both the den- 
sity of resonant absorption and the amplitude of the X mode opacity 
resonance induced by vacuum corrections decline with photon en- 
ergy. Radiation of sufficiently low frequency will have only modest 
interaction with the opaque layer and will retain the plasma dom- 
inated photosphere (r^, p <C 1). This may be contrasted with ra- 
diation which decouples from within the opaque layer because of 
the large optical depth accumulated at p res . For the 10 14 G models 
illustrated, this transition occurs at ~ 1.5 — 3.0 keV. Figure 1 1 2i 
shows the decoupling density of these strongly magnetized mod- 
els. At energies > 3 keV, the X mode radiation decouples from 
somewhat cooler strata than would be found were VP effects ne- 
glected; no radiation decouples from the relatively tenuous vacuum 
dominated plasma (t^P 3> 1) and the response properties charac- 
teristic of this regime have no influence on either the polarization 
or angular distribution of the spectra (see i|2.2> . 

Thermal absorption in the opaque heating layer has an ancil- 
lary effect on the spectrum at long wavelengths which decouple 
below the layer, where the flux amplitude is enhanced by the re- 
radiated emission. The formation depth of the proton cyclotron line 
becomes comparable to the depth of the heating layer, which re- 
duces the effective width of the proton cyclotron line in the 10 14 G 
models. For models log T e g > 6.5 the cyclotron absorption line 
is preserved but narrower than expected from the equivalent model 
without vacuum corrections (Fig. I13i . For cooler T e ff, the forma- 
tion depth of the line coincides with larger T than for neighboring 
frequencies, leading to an apparent emission feature. Regardless of 
T e ff, the reduced effective width of the line results from enhanced 
thermal emission from the opaque layer. That this tranformation 
occurs over a range of just 2.5 T e g suggests that this feature may 
be of diagnostic value in thermal spectra of some anomalous X-ray 
pulsars. Figure i 14i summarizes the evolution in cyclotron temper- 
ature with T e ff . 



6 DISCUSSION 

Observations of thermal spectra from cooling isolated NSs provide 
a direct probe of their surface properties, assuming these data can 
be adequately interpreted with a realistic spectral model. The sur- 
face magnetic field strength can be measured through the detec- 
tion of cyclotron features. Strong atomic transitions constrain the 
ratio M/R, while an additional measure of the gravitational poten- 
tial (e.g. through line profiles) could restrict M and R individually 
(Paerels 1997); line broadening mechanisms and variation of sur- 
face properties pose additional difficulty for interpretation of these 
properties. With the exception of IE 1207.4-5209 for which the na- 
ture of absorption features is presently uncertain, the failure to mea- 
sure any spectral features in the Chandra and XMM observations 
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of thermally emitting NS is notorious. Non-detection of features 
in NS surface emission tends to favor light element light element 
atmosphere models; the Fe atmospheres of Raia eopal et alJ J 19971) 
for example predict strong absorption edges at soft X-ray energies 
which should be unambiguous in the data. The absence of cyclotron 
features may exclude ranges of surface field strength but, for well 
magnetized stars, such features can be smoothed by realistic dis- 
tributions of the surface temperature and field iZane et al. 2001). 
Measurement of the integrated net polarization at X-ray energies 
for cooling, isolated NS may provide enough information to deduce 
the surface magnetic field strength where spectra l features remain 
unresolved or are dest royed by other processes l Pavlov & Zavlin 
bOOOHHevl et all2003l) . 

We have presented a flexible, efficient and robust method for 
computation of light element thermal NS spectra and atmospheric 
structures. The radiative processes in these models are the simplest 
required to describe the X-ray properties of cooling light element 
NS atmospheres for a broad range of stellar parameters. A number 
of refinements to the model physics and distinguishing elements of 
the present work are described below. 

Ionization equilibrium — In strong magnetic fields the ground 
state binding energy of stationary Hydrogen atoms grows asymp- 
totically with ~ In 2 B c Ryd where B c = B/2.35 x 10 9 G iLouded 
1959), and bound state transitions may become important to the 
X-ray spectrum of NS at 10 12 G. Motional effects of the atoms 
within strong magnetic fields leads to the formation of decentered 
states , and yields substan tially altered cross sections (Potekhin 
1998). IPavlov et alJ ^995) have shown that even a modest neu- 
tral plasma component can alter the spectrum owing to bound- 
free transitions from these tightly bound states. Moreover, pres- 
sure destruction is less effective in strong magnetic fields owing 
to the reduce phase space volume accessible to the electron distri- 
bution. An occupation probability (OP) formalism has been studied 
in detail for both magnetized Potekhin et alJ 1 19991) and unmagne- 
tized Hydrogen plasma I Potekhinl l 19961) . This free-energy picture 
is parametrized by the binding energies and mean atomic sizes, and 
provides a natural equation of state governing the pressure destruc- 
tion of weakly bound states and truncation of the partition function. 
Steady progress has been make in the EOS for magnetized hydro- 
gen, and a recent report describs an effort to integrate the ionization 
equilibrium of Potekhin & Chabrier 1 2002) with the Lucy-Unsold 
temperature correction scheme (Ho et al. 2002). 

It is more challenging to incorporate non-LTE processes to the 
statistical equilibrium, by which the bound state distributions are 
determined through detailed balance of the collisional and radiative 
rate coefficients. The radiative rates can be evaluated by straightfor- 
ward integration given (e.g.) the bound-free cross sections and the 
state of the radiation field. The collisional transition rates must also 
be specified, unlike LTE in which it is sufficient to consider ra- 
tios of these rates and which are given by Boltzmann factors; these 
rates, and preferably their derivatives, must be computed or other- 
wise tabulated for the spectrum of bound states and continuum. The 
CL technique is able to respond to neutral or incompletely ionized 
component better than Lucy-Unsold upon simultaneous integration 
of the statistical equilibrium and radiative transfer, but may require 
adaptive depth points around the transition zones between succes- 
sive ionization states or near line formation depths. Much simpler 
would be inclusion of OP AL opac ities in the present code (see e.g. 
IWerner & Deeded feOOOl) or lGansicke et alJ 120021) ) but these tab- 
ulations are not appropriate for strongly magnetized atmospheres. 
The ionized model calculations described in this paper underutilize 
the advantages afforded by the CL method. 



High density effects — Two phenomena related to dense 
plasma are not accounted for in the opacity formulation of [J2| Col- 
lective plasma effects are encountered where v e > 1 strongly in- 
fluence the optical and near UV flux of thermal models with ra- 
dio pulsar type fields (Fig. [6j; the effects extend to X-ray ener- 
gies for the magnetar models (Fig. ll2t . In this regime, the collec- 
tive response of the plasma should strongly attenuate propagating 
waves, and the assumption of radiative transport may be incorrect 
to long wavelength radiation. This is plain in the magnetar spectra, 
but is also relevant at optical wavelengths for very cool stars (e.g. 
RXJ 1856.35) and aging RPs like Geminga. Collective phenom- 
ena may be included by way of a structure factor in the scattering 
opacity as a first step toward evaluating corrections to the equilib- 
rium thermal structure. Secondly, plasma is tightly coupled when 
T = (47tn e /3) 1/3 e 2 /(kT) > 1, but the Debye-Huckel approx- 
imation 1481 is only a perturbative correction. A more plausible 
model for Pcoui could be implemented, although errors in the cou- 
pling occur at the very deepest model strata of the most strongly 
magnetized atmospheres where an improved calculation is likely to 
correct only the X mode intensities in the spectral tail. 

Ionic cyclotron resonances — Our opacity calculation is de- 
rived from the self-consistent inclusion of ions in the plasma di- 
electric J2j; energy dependencies which scale with magnetic field 
enter the opacity through the eigenvalues of the polarization ten- 
sor 1201 . Some authors favor a calculation of these eigenvalues in 
which electrons and protons are distinct plasma constituents, with 
effective eigenvalues (c.f. lZane etai]j2000]) : lHo & LaiH200ll) ): 

J 2 1 , (m e /m p ) 2 

As discussed in fl2.ll the principal source of ion cyclotron ab- 
sorption resides in resonant behavior of the polarization vector it- 
self and, if the vectors e J p are computed consistently as described 
in the discrepancy of no particular consequence for radiation 
uj > lj c ,p', the opacity to longer wavelength radiation however ac- 
quires a scaling factor oc from eqn <20k which is steeper than 
the ~ u^ 1 scaling of <6H . Results calculated in the approximation 
<6H are potentially unreliable at X-ray energies for B > 10 14 G, 
and the discrepancy is compounded for admixtures of helium. 

Warm plasma and other opacity effects — The opacity model 
of |2|is derived assuming "cold" electron plasma plus vacuum and 
neglects plasma disp ersion which is o f interest near the electron 
cyclotron line ( Kirk & Meszaros 1980). We have also neglected 
Comptonization which affe cts the am plitude of the high energy 
spectral cutoff iKaminker et"ai]|l983l) . For the NS model temper- 
atures discussed herein, a fully relativistic treatment of the plasma 
response is unnecessary (Meszaros 1992). Calculation of the mag- 
netic Gaunt factors at transitional field strengths I ~ 1 are incom- 
plete as described in <j2]and could be improved by using the calcu- 
lation of lPavlov & Panovl ll976l) . 
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Figure 1. Total emergent flux profiles for unmagnetized Hydrogen plasma 
atmosphere (solid curves), and weakly magnetized plasma log B = 
8.0(1.0)10.0 G, at three effective temperatures log T eff = 5.6(0.4)6.4. 
For logB = 10.0 the cyclotron absorption encroaches on the soft X-ray 
spectrum, and at lower T e g becomes comparable to the peak emission, 
yielding enhanced optical and near UV amplitudes. With the exception of 
the coolest logf? = 10.0 atmospheres, the model spectra satisfy a single 
displacement rule E™ ax ~ 4.7kT eff . 
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Figure 2. Spectral envelopes for log B = 10.0 at three flux temperatures, 
log T e g = 5.6(0.4)6.4; resonant absorption significantly affects the spec- 
trum where the electron cyclotron feature is nearly coincident with peak 
emission energy. Only X-mode radiation is sensitive to electron cyclotron 
absorption, and energy thus absorbed is reemitted at longer wavelenths. 



Figure 4. Temperature profiles for several converged, weakly magnetized 
models, illustrating the variation in asymptotic surface temperature which 
results from energy redistribution from cyclotron absorption in a radiative 
equilibrium atmosphere. Curves are labelled by B = (unmagnetized) 
or log B. At finite B, departures from the unmagnetized thermal structure 
occur at deeper strata for increasing field strength. 
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Figure 3. Polarization fraction (upper), photospheric temperature (middle) 
and decoupling density (lower panel) for the model spectra of FigEI Reso- 
nant absorption plays a more significant role at cooler temperatures, where 
the peak thermal energy approximately coincides with the cyclotron line. 
Thomson scattering of high energy photons is energy independent; lower 
energy photons which are predominantly absorbed decouple from similar 
densities regardless of T e g . 
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Figure 5. Model spectra for logB = 12.0 and log T eff = 5.6(0.4)6.4 
including the partial polarized flux contributions. The broad electron cy- 
clotron line is centered off diagram at E ~ 11.6 kcV, but its influence is 
seen in the relative amplitude of emission in the two modes. 
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Figure 6. Polarization and decoupling density for the model spectra of 
Fig J5j- Radiation near the electron feature have large negative H(E) ow- 
ing to the relatively large O mode amplitude (upper panel). The proton line 
demonstrates modest variation in net polarization. The v~ 1 dependence of 
the X mode opacity is more apparant than for the 10 10 G models (lower 
panel). Models with higher T e g decouple at greater p. Measurement of 
Tl(E) in two bands may be sufficient to diagnose T e g for surface fields 
B ~ 10 11 — 10 13 G. 
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Figure 7. Equilibrium temperature profiles on the mean Rosseland depth 
scale for log T e g = 5.7, 6.4 with intermediate magnetic fields log B = 
11.0(1.0)13.0. As B increases, a greater fraction of the total flux is emitted 
in the (transparant) X mode from deep, hot strata. The asymptotic surface 
temperature declines with increasing field strength owing to the depleted 
intensity in the more opaque O mode. 



Figure 8. Angular intensity profiles near the spectral peak for log T e jj = 
6.0, log B = 12.0, and four choices of magnetic field orientation. For each 
curve the vectors B,k,n are coplanar (rp = 0). Intensity along the field 
vector declines with increasing inclination by decoupling from lower tem- 
peratures in higher strata. The pencil beam has approximately equal contri- 
butions from both modes while the fan emission is highly polarized. 
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Angular intensity profiles for B = 


10 12 , T eff 
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b = 30° at a photon energy near the spectral peak for six choices of model 
azimuths. The pencil amplitude declines and acquires finite polarization for 
<t> > and vanishes for <f> > 9 V and the fan component dominates for large 

4>. 
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Figure 10. Total flux profiles for log B = 12.0, log T eff = 5.6(0.4)6.4 
and inclinations 6 = 0, 30, 60, 90. Spectra for T e g > 6.0 are softer at 
high energies for larger field inclinations. The converged T(m) for these 
models do not reveal large variation with b for a given T e g. 
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Figure 12. The decoupling density for B = 10 14 G at log T eff = 6.7, 
illustrating the transition to photosphere formation within the vacuum reso- 
nance layer above 1.5 keV. The X-mode flux below 0.25 keV is dominated 
by collective plasma effects. 
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Figure 11. Thermal structure of the magnetar models B = 10 G for a 
range of effective temperatures log T e g = 6.2(0.1)6.7. The effect of the 
vacuum resonance on the transparant X mode is clearly seen in the form 
of the heating layer. Small "jumps" in the temperature profiles arise from 
abrupt variation in opacity at the adaptive depth points, particularly for en- 
ergies E > 1 keV. The dashed profile is that of the converged model for 
pure plasma at log T eff = 6.7 and B = 10 14 G. 



Figure 13. Total flux profiles for log T eff = 6.3(0.2)6.7 at 10 14 G, com- 
paring the vacuum polarized models (solid curves) to their ordinary plasma 
equivalents (long dash). Reemission of thermal radiation absorbed in the 
opaque layer (vacuum resonance) effectively "fills in" the proton cyclotron 
line in hot models, and inverts the line profile in relatively cool models. 
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Figure 14. The flux profiles for the magnetar models B = 10 G for the 
effective temperatures log T e g = 6.2(0.1)6.7. The width of the proton 
cyclotron line is reduced by the incidental effects of vacuum polarization 
on the atmospheric structure (Fig ll II . This phenomenon is ultimately re- 
sponsible for line "inversion" in cooler models. 



